E013 Pho4 chimera activity analysis, new host

require(plotly)

Loading required package: plotly
Loading required package: ggplot2
Keep up to date with changes at https://www.tidyverse.org/blog/
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

require(tidyverse)

Loading required package: tidyverse
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ tibble  3.1.6     ✓ dplyr   1.0.7
✓ tidyr   1.1.4     ✓ stringr 1.4.0
✓ readr   2.1.1     ✓ forcats 0.5.1
✓ purrr   0.3.4     
── Conflicts ───────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks plotly::filter(), stats::filter()
x dplyr::lag()    masks stats::lag()

require(cowplot)

Loading required package: cowplot

Goal

• Analyze the latest batches of flow cytometry data to determine the contribution of different (matching) parts of ScPho4 and CgPho4 to their difference in Pho2 dependence.
• Develop an analysis pipeline to perform QC, correction (if needed) and plotting the results.

Data

Merge the three batches

raw <- list(
  "03/15/22" = read_tsv("../data/20220315-EO-chimera-batch-1/20220329-gated-median-out.txt", show_col_types = FALSE),
  "03/22/22" = read_tsv("../data/20220322-EO-chimera-batch-2/20220327-gated-median-out.txt", show_col_types = FALSE),
  "03/29/22" = read_tsv("../data/20220329-EO-chimera-batch-3/20220330-gated-median-out.txt", show_col_types = FALSE)
)
dat <- bind_rows(raw, .id = "date") %>% 
  #separate(col = well, sep = 1, into = c("row", "col")) %>% 
  separate(col = sample, sep = "-", into = c("plasmid", "host", "rep"))

QC

1. Number of events per sample

dat %>% ggplot(aes(x = n_induction)) + geom_histogram(bins = 30, fill = "forestgreen") + facet_grid(host~date) + theme_gray(base_size = 14) #+ panel_border()

Discussion

1. Majority of the experimental wells (not blank) have > 10000 events, with the exception of a few samples
2. Remove the blank wells and proceed with QC

2. Background subtraction

Check the background fluorescence levels across different days.

dat %>% filter(host == "156") %>% 
  pivot_longer(BL1.H:nRFP, names_to = "parameter", values_to = "intensity") %>% 
  mutate(parameter = ordered(parameter, levels = c("BL1.H", "YL2.H", "nGFP", "nRFP"))) %>% 
  ggplot(aes(x = date, y = intensity)) + geom_point(aes(shape = well)) +
  stat_summary(geom = "pointrange", color = "red", position = position_nudge(x = 0.1)) +
  facet_wrap(~parameter, scale = "free_y") +
  theme_bw(base_size = 14)

No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`

03/29/22 D9 showed higher background fluorescence than the others. after checking the gating result, I found that NA-156-2 sample in D9 had a slightly wider distribution, and the 2d cluster gate is also wider than the other two samples. The cell size normalized intensity doesn’t show the same pattern, suggesting that at least some of the difference could be due to cell size variation. Try just averaging the three blank wells and perform the background removal.

Subtract the background

bg <- dat %>% 
  filter(host == "156") %>% 
  group_by(date) %>% 
  summarize(across(FSC.H:nRFP, ~ round(mean(.x),1))) %>% 
  column_to_rownames(var = "date")

dat1 <- dat %>% 
  filter(host != "blank") %>% 
  select(date:host, FSC.H:nRFP) %>%
  mutate(
    BL1.H = BL1.H - bg[date, "BL1.H"],
    YL2.H = YL2.H - bg[date, "YL2.H"],
    nGFP = nGFP - bg[date, "nGFP"],
    nRFP = nRFP - bg[date, "nRFP"],
  )

Double check that the subtraction worked correctly

dat1 %>% filter(host == "156") %>% 
  pivot_longer(BL1.H:nRFP, names_to = "parameter", values_to = "intensity") %>% 
  mutate(parameter = ordered(parameter, levels = c("BL1.H", "YL2.H", "nGFP", "nRFP"))) %>% 
  ggplot(aes(x = date, y = intensity)) + geom_point(aes(shape = well)) +
  stat_summary(geom = "pointrange", color = "red", position = position_nudge(x = 0.1)) +
  facet_wrap(~parameter, scale = "free_y") +
  theme_bw(base_size = 14)

No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`

# remove the yH156 samples
dat1 <- dat1 %>% 
  filter(host != "156") %>% 
  mutate(host = ordered(host, levels = c("555", "373"), labels = c("PHO2", "pho2∆")))

Export the background-subtracted values for later use with the shinyapp

write_tsv(dat1, "../shinyapp/20220412-bg-sub-data.tsv")

3. Consistency across plates.

Check the background-subtracted fluorescence values for the host strains (yH373 and yH555) as well as the positive control strain (pH188 in yH373 or yH555 backgrounds), both of which are present on each plate.

dat1 %>% filter(plasmid %in% c("188", "NA")) %>% 
  pivot_longer(BL1.H:YL2.H, names_to = "parameter", values_to = "intensity") %>% 
  #mutate(parameter = ordered(parameter, levels = c("BL1.H", "YL2.H", "nGFP", "nRFP"))) %>% 
  ggplot(aes(x = plasmid, y = intensity, color = date)) + 
  geom_point(position = position_dodge(0.5)) + 
  scale_color_brewer(type = "qual", palette = 2) +
  facet_grid(parameter~host, scale = "free_y") +
  theme_bw(base_size = 14) +
  theme(axis.text.x = element_text(face = 3))

1. The host strain fluorescence levels appear to be consistent across the three days. Additionally, their GFP levels are zero, as expected.
2. The positive control strain has strong BL1.H and correspondingly has strong YL2.H.

The size normalized values show slightly lower variance:

dat1 %>% filter(plasmid %in% c("188", "NA")) %>% 
  pivot_longer(nGFP:nRFP, names_to = "parameter", values_to = "intensity") %>% 
  #mutate(parameter = ordered(parameter, levels = c("BL1.H", "YL2.H", "nGFP", "nRFP"))) %>% 
  ggplot(aes(x = plasmid, y = intensity, color = date)) + 
  geom_point(position = position_dodge(0.5)) + 
  scale_color_brewer(type = "qual", palette = 2) +
  facet_grid(parameter~host, scale = "free_y") +
  theme_bw(base_size = 14) +
  theme(axis.text.x = element_text(face = 3))

4. Well position

Rationale

• In my original plate design, I placed a positive control strain (envisioned CgPho4-mNeon in PHO2 background for example) in the control columns (1, 5, 9) on every other row (A, C, E, G). The reason for this is to use it to identify any well position effect on the fluorescence readings.
• In Emily’s implementation, she instead put a pair of strains, namely pH188-yH323 and pH188-yH555 in these wells. This means I cannot use the positive control wells exactly the way as I designed them. But I can still use them to spot any trend.

dat1 %>% 
  filter(plasmid == "188") %>% 
  separate(well, into = c("row", "col"), sep = 1) %>% 
  pivot_longer(BL1.H:YL2.H, names_to = "parameter", values_to = "intensity") %>% 
  ggplot(aes(x = row, y = intensity, group = date)) + 
  geom_point(aes(color = date), size = 1) + 
  #geom_line(aes(color = date), size = 0.4) +
  stat_summary(aes(color = date), geom = "line", fun = mean) +
  scale_color_brewer(type = "qual", palette = 2) +
  facet_grid(parameter~host, scales = "free_y") +
  theme_bw(base_size = 14)

• no obivous trend between the rows (or columns, not shown)
• a clear correlation between YL2.H and BL1.H. at least part of this is due to the cell size differences – see size normalized data below:

dat1 %>% 
  filter(plasmid == "188") %>% 
  separate(well, into = c("row", "col"), sep = 1) %>% 
  pivot_longer(c(YL2.H, nRFP), names_to = "parameter", values_to = "intensity") %>% 
  ggplot(aes(x = row, y = intensity, group = date)) + 
  geom_point(aes(color = date), size = 1) + 
  #geom_line(aes(color = date), size = 0.4) +
  stat_summary(aes(color = date), geom = "line", fun = mean) +
  scale_color_brewer(type = "qual", palette = 2) +
  facet_grid(parameter~host, scales = "free_y") +
  theme_bw(base_size = 14)

• What we care about is the ratio between Pho4-mNeon and PHO5pr-mCherry

dat1 %>% 
  filter(plasmid == "188") %>% 
  separate(well, into = c("row", "col"), sep = 1) %>% 
  ggplot(aes(x = BL1.H, y = YL2.H)) + 
  geom_point(aes(shape = date, color = row), size = 3) + 
  stat_smooth(method = "lm", se = FALSE, color = "gray50", size = 0.5) +
  scale_color_viridis_d() +
  scale_shape_manual(values = 15:17) +
  #scale_x_log10() + scale_y_log10() +
  #facet_grid(parameter~host, scales = "free_y") +
  theme_bw(base_size = 14)

`geom_smooth()` using formula 'y ~ x'

• There is variation in BL1.H and correspondingly in YL2.H, but the ratio between the two are very consistent

Let’s check both ScPho4 (pH194) and CgPho4 (pH188) to see if their behaviors are consistent across the days.

dat1 %>% 
  filter(plasmid %in% c("188", "194")) %>% 
  mutate(`YL2.H/BL1.H` = YL2.H/BL1.H,
  #mutate(`YL2.H/BL1.H` = nRFP/nGFP,
         Pho4 = factor(plasmid, levels = c("194", "188"), 
                       labels = c("ScPho4", "CgPho4"))) %>% 
  separate(well, into = c("row", "col"), sep = 1) %>% 
  ggplot(aes(x = host, y = `YL2.H/BL1.H`)) + 
  geom_point(color = alpha("gray50", 0.9), position = position_jitter(0.2)) + 
  stat_summary(fun.data = "mean_se", geom = "pointrange", color = "red") +
  facet_grid(Pho4~date) +
  theme_bw(base_size = 14)

5. High variance samples

Summarize the background subtracted data by calculating the means and cv for each strain.

dat2 <- dat1 %>% 
  select(-nGFP, -nRFP) %>% 
  pivot_longer(FSC.H:YL2.H, names_to = "parameter", values_to = "intensity") %>% 
  group_by(date, plasmid, host, parameter) %>% 
  summarize(
    n = n(),
    mean = num(mean(intensity), digits = 0),
    cv = num(sd(intensity)/mean(intensity), digits = 2)
  ) %>% 
  arrange(desc(cv))

`summarise()` has grouped output by 'date', 'plasmid', 'host'. You can override using the `.groups` argument.

dat2 %>% 
  filter(plasmid != "NA", parameter != "FSC.H") %>% 
  ggplot(aes(x = cv)) + geom_histogram(aes(y = ..density../20)) + stat_ecdf() + 
  geom_hline(yintercept = 0.8, linetype = 2) +
  ylab("cumulative density") + 
  facet_wrap(~parameter) +
  theme_cowplot()

Don't know how to automatically pick scale for object of type pillar_num/pillar_vctr/vctrs_vctr/double. Defaulting to continuous.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

The histogram’s y-axis is not shown. The line graph represents the empirical CDF, and the dotted horizontal line is at 80%. GFP is more variable than RFP, likely because the absolute values of the former is lower. For both, ~80% of the samples have a CV < 20%.

Do the same for the cell size-normalized values

ndat2 <- dat1 %>% 
  select(-BL1.H, -YL2.H) %>% 
  pivot_longer(nGFP:nRFP, names_to = "parameter", values_to = "intensity") %>% 
  group_by(date, plasmid, host, parameter) %>% 
  summarize(
    n = n(),
    mean = num(mean(intensity), digits = 0),
    cv = num(sd(intensity)/mean(intensity), digits = 2)
  ) %>% 
  arrange(desc(cv))

`summarise()` has grouped output by 'date', 'plasmid', 'host'. You can override using the `.groups` argument.

ndat2 %>% 
  filter(plasmid != "NA", parameter != "FSC.H") %>% 
  ggplot(aes(x = cv)) + geom_histogram(aes(y = ..density../20)) + stat_ecdf() + 
  geom_hline(yintercept = 0.8, linetype = 2) +
  ylab("cumulative density") + 
  facet_wrap(~parameter) +
  theme_cowplot()

Don't know how to automatically pick scale for object of type pillar_num/pillar_vctr/vctrs_vctr/double. Defaulting to continuous.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

CV is smaller on the cell size normalized fluorescence values.

Main

Plotting functions

Here I’d like to develop a series of plotting functions that take the names or any part of the Pho4 chimera annotation as input and plot their results in a variety of ways Load the Pho4 plasmid information and merge with the reshaped data

meta <- read_tsv("../data/20220330-chimera-Pho4-makeup.tsv", col_types = "cccc")

Plotting function

# subset data
tmp <- meta %>% 
  filter(symbol %in% c(refs,test)) %>% 
  inner_join(dat1, by = "plasmid") %>% 
  mutate(symbol = factor(symbol, levels = c(refs, test)))

Error: Problem with `mutate()` column `symbol`.
ℹ `symbol = factor(symbol, levels = c(refs, test))`.
x factor level [3] is duplicated
Run `rlang::last_error()` to see where the error occurred.

Plot individual components

tmp %>% 
  pivot_longer(nGFP:nRFP, names_to = "parameter", values_to = "intensity") %>% 
  mutate(parameter = ordered(parameter, levels = c("nRFP", "nGFP"))) %>% 
  #pivot_longer(BL1.H:YL2.H, names_to = "parameter", values_to = "intensity") %>% 
  #mutate(parameter = ordered(parameter, levels = c("YL2.H", "BL1.H"))) %>% 
  ggplot(aes(x = symbol, y = intensity, group = host)) + 
  geom_bar(aes(fill = host), width = 0.5, alpha = 0.8,
           stat = "summary", fun = "mean", position = position_dodge(0.5)) +
  geom_point(aes(color = host), size = 1,
             position = position_jitterdodge(dodge.width = 0.5)) + 
  scale_color_manual(values = c("PHO2" = "gray20", "pho2∆" = "gray40")) +
  scale_fill_viridis_d(begin = 0.2, end = 0.6) +
  #stat_summary(fun = "mean", geom = "crossbar", color = "red", width = 0.25,
  #             position = position_dodge(0.75), ) +
  facet_wrap(~parameter, scale = "free_y", ncol = 1) +
  xlab("Pho4 chimera") + expand_limits(y = 0) +
  theme_gray(base_size = 14) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, family = "mono"))

tmp %>% 
  #mutate(`YL2.H/BL1.H` = YL2.H/BL1.H,
  mutate(`R/G` = nRFP/nGFP) %>% 
  ggplot(aes(x = symbol, y = `R/G`, group = host)) + 
  geom_bar(aes(fill = host), width = 0.5, alpha = 0.8,
           stat = "summary", fun = "mean", position = position_dodge(0.5)) +
  geom_point(aes(color = host), position = position_jitterdodge(dodge.width = 0.5)) + 
  scale_color_manual(values = c("PHO2" = "gray20", "pho2∆" = "gray40")) +
  scale_fill_viridis_d(begin = 0.2, end = 0.6) +
  #stat_summary(fun = "mean", color = "red", geom = "crossbar", width = 0.2,
  #             position = position_dodge(0.75), ) +
  theme_gray(base_size = 14) + xlab("Pho4 chimera") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, family = "mono"))

Next idea is to plot the ratio between the mean R/G for PHO2 and mean R/G for pho2∆, over the mean R/G for PHO2

First summarize the data.

datsum <- dat1 %>%
  filter(plasmid != "NA") %>% 
  mutate(`R/G` = YL2.H/BL1.H, `nR/G` = nRFP/nGFP) %>% 
  group_by(date, plasmid, host) %>% 
  summarize(across(c(BL1.H, nGFP, YL2.H, nRFP, `R/G`, `nR/G`), mean), .groups = "drop") %>% 
  pivot_wider(names_from = host, values_from = BL1.H:`nR/G`) %>% 
  mutate(`pho2∆/PHO2` = `R/G_pho2∆`/`R/G_PHO2`,
         `n.pho2∆/PHO2` = `nR/G_pho2∆`/`nR/G_PHO2`)

# useful to set a flag for low activity mutants
low.act.th <- datsum %>% filter(plasmid == "194") %>% summarize(across(.cols = where(is.numeric), .fns = mean)) %>% unlist()

Then extract the subset for plotting

tmpsum <- meta %>% 
  filter(symbol %in% c(refs,test)) %>% 
  inner_join(datsum, by = "plasmid") %>% 
  mutate(symbol = factor(symbol, levels = c(refs, test)))

Design the plot

tmpsum %>% 
  mutate(Activity = ifelse(`R/G_PHO2` < 2*low.act.th["R/G_pho2∆"], "low", "pass")) %>% 
  ggplot(aes(x = symbol, y = `pho2∆/PHO2`)) + 
  geom_col(aes(group = date, fill = `R/G_PHO2`), width = 0.75, color = "gray50",
           position = position_dodge(0.9)) +
  #geom_point(aes(color = host), position = position_jitterdodge(dodge.width = 0.5)) + 
  scale_fill_gradient2("Activity") +
  #scale_color_manual(values = c(alpha("black",0), "red3")) +
  #stat_summary(fun = "mean", color = "red", geom = "crossbar", width = 0.2,
  #             position = position_dodge(0.75), ) +
  facet_grid(.~Activity, scales = "free_x", space = "free_x", labeller = "label_both") +
  theme_bw(base_size = 14) + xlab("Pho4 chimera") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, family = "mono"))

X-Y plot

p3 <- tmpsum %>% 
  mutate(`nR/G_PHO2` = signif(`nR/G_PHO2`, digits = 2),
         `nR/G_pho2∆` = signif(`nR/G_pho2∆`, digits = 2)) %>% 
  ggplot(aes(x = `nR/G_PHO2`, y = `nR/G_pho2∆`, label = symbol)) + 
  geom_point(size = 2.5) + geom_abline(slope = 1) +
  theme_gray(base_size = 14)
ggplotly(p3, tooltip = c("label", "x", "y"))